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Abstract 

This paper describes the implementation of opti- 
mization techniques based on control theory for air- 
foil design. In previous studies (6, 7] it was shown 
that control theory could be used to devise an ef- 
fective optimization procedure for two-dimensional 
profiles in which the shape is determined by a con- 
formal transformation from a unit circle, and the 
control is the mapping function. The goal of our 
present work is to develop a method which does not 
depend on conformal mapping, so that it can be ex- 
tended to treat three-dimensional problems. There- 
fore, we have developed a method which can address 
arbitrary geometric shapes through the use of a fi- 
nite volume method to discretize the potential flow 
equation. Here the control law serves to provide 
computationally inexpensive gradient information to 
a standard numerical optimization method. Results 
are presented, where both target speed distributions 
and minimum drag are used as objective functions. 

Nomenclature 

Aij grid transformation coefficients 
b design variable 
B generic co-state variable 
c speed of sound 

C bounding surface of flowfield domain on airfoil 
Cd coefficient of drag 
Ct coefficient of lift 
C p coefficient of pressure 
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Cp coefficient of pressure for sonic flow 
d modulus of conformal mapping transformation 
D flowfield domain 
/ control law, control parameter 
G generic flowfield variable 
h Jacobian of generalized transformation 
H grid transformation matrix 
/ cost function 
J grid transformation matrix 
m number of flowfield evaluations per line search 
M local Mach number 
Afoo Mach number at infinity 
n number of design variables 
p pressure 

P grid perturbation variational 
q speed 

qd desired speed 
speed at infinity 

Q grid perturbation variational 
R generic governing equation for flowfield 
3 arc length along airfoil surface 
5 modified boundary condition for 
<i , t-i exponents on basis functions 
u, v Cartesian velocity components 
U, V contra-variant velocity components 
x, y Cartesian coordinates 
X physical position in flowfield 
rj body fitted coordinates 
a angle of attack 
<j> velocity potential 
7 ratio of specific heats 

co-state variable for potential flow 
p density 

9 angle around circle 
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Formulation of the design 
problem as a control problem 

Ultimately, the designer seeks to optimize the geo- 
metric shape of a configuration taking into account 
the trade-offs between aerodynamic performance, 
structure weight, and the requirement for internal 
volume to contain fuel and payload. The subtlety 
and complexity of fluid flow is such that it is un- 
likely that repeated trials in an interactive analy- 
sis and design procedure can lead to a truly opti- 
mum design. Progress toward automatic design has 
been restricted by the extreme computing costs that 
might be incurred from brute force numerical opti- 
mization. However, useful design methods have been 
devised for various simplified cases, such as two- 
dimensional airfoils in viscous flows [13] and wings 
in invicid flows. The computational costs for these 
methods result directly from the vast number of flow 
solutions that are required to obtain a converged de- 
sign. 

Alternatively, it has been recognized that the de- 
signer generally has an idea of the kind of pres- 
sure distribution that will lead to the desired per- 
formance. Thus, it is useful to consider the inverse 
problem of calculating the shape that will lead to a 
given pressure distribution. The method is advan- 
tageous, since only one flow solution is required to 
obtain the desired design. Unfortunately, a physi- 
cally realizable shape may not necessarily exist, un- 
less the pressure distribution satisfies certain con- 
straints, thus the problem must be very carefully 
formulated. 

The problem of designing a two-dimensional pro- 
file to attain a desired pressure distribution was first 
studied by Lighthill, who solved it for the case of in- 
compressible flow with a conformal mapping of the 
profile to a unit circle [9]. The speed over the profile 
is 

« = - d |Vrf| 

where <f> is the potential which is known for incom- 
pressible flow and d is the modulus of the mapping 
function. The surface value of d can be obtained 
by setting q = qj, where qj is the desired speed, and 
since the mapping function is analytic, it is uniquely 
determined by the value of d on the boundary. A so- 


lution exists for a given speed qoo at infinity only if 



and there are additional constraints on q if the pro- 
file is required to be closed. 

The difficulty that the objective may be unattain- 
able can be circumvented by regarding the design 
problem as a control problem in which the control 
is the shape of the boundary. A variety of alterna- 
tive formulations of the design problem can then be 
treated systematically within the framework of the 
mathematical theory for control of systems governed 
by partial differential equations [10]. This approach 
to optimal aerodynamic design was introduced by 
Jameson [6, 7], who examined the design problem 
for compressible flow with shock waves, and devised 
adjoint equations to determine the gradient for both 
potential flow and also flows governed by the Euler 
equations. More recently Ta’asan, Kuruvila, and 
Salas, implemented a one shot approach in which 
the constraint represented by the flow equations is 
only required to be satisfied by the final converged 
solution [16]. Pironneau has also studied the use of 
control theory for optimum shape design of systems 
governed by elliptic equations [12]. 

Suppose that the boundary is defined by a func- 
tion /(x), where x is the position vector, and the 
desired objective is measured by a cost function I. 
This may, for example, measure the deviation from 
a desired surface pressure distribution, but it can 
also represent other measures of performance such 
as lift and drag. Thus, the problem is recast into a 
numerical optimization procedure in which the com- 
putationally expensive finite difference gradient of 
the objective function with respect to the control or 
design variables is replaced by the first variation in 
the cost function. Suppose that a variation 5/ in 
the control produces a variation SI in the cost. 

For flow about an airfoil the aerodynamic proper- 
ties which define the cost function are functions of 
the flow-field variables ( G ) and the physical location 
(*); 

I = I(G,X). 

Thus, 

s, = ld ,a+ § sx - (1 > 

As pointed out by Baysal and Elesbaky [1] each term 
in (1), except for 5G, can be easily obtained. 
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and jx can be obtained directly without a flowfield 
evaluation since they are partial derivatives. 8X can 
be determined for each design variable by succesive 
grid generation so long as this cost is significantly 
less then the cost of the flow solution. 6G tradition- 
ally requires a flowfield evaluation for each design 
variable independantly (brute force gradient meth- 
ods). Here we introduce the governing equations of 
the flowfield as a constraint in such a way that the 
flowfield evaluation is eliminated from the expression 
for the gradient. Like the cost function, the govern- 
ing equation R is a function of G and X within the 
flowfield domain D , 

R = R{G,X) = 0. 


Thus, 


SR ‘^ SG+ § SX = 0 - 


( 2 ) 


dl 8R 

dG 3 dG 

= 0 , 

( 3 ) 

§c 6x ~ B 1 

{.dx 6X ■ 

( 4 ) 


Next, using a Lagrange Multiplier B we have 

51 = %d SG + JZ SX ~ B + JZ 6X 

Choosing B such that, 


gives 


The advantage is that (4) is independent of 5G , 
therefore additional flowfield evaluations are not re- 
quired to determine the gradient of I with respect 
to any number of design variables. The main cost is 
in solving the adjoint equation (3). In general, this 
adjoint problem is about as complex as a flow solu- 
tion. If the number of design variables is large then 
the cost differential between one adjoint solution and 
many flowfield evaluations required to determine the 
gradient becomes compelling. Instead of intorduc- 
ing a Lagrange multiplier, B, one can solve (2) for 
8G as 

'dR' 
dG . 

and insert the result in (1). This is the implicit 
gradient approach, which is essentially equivalent to 
the control theory approach as has been pointed out 
by Shubin and Frank [14, 15]. 

After making such a modification, the gradient 
can be recalculated and the process repeated to fol- 
low a path of steepest descent until a minimum is 


8G = - 


i-i 


™SX 

dx 6X ' 


reached. In order to avoid violating constraints, such 
as a minimum acceptable wing thickness, the gradi- 
ent may be projected into the allowable subspace 
within which the constraints are satisfied. In this 
way one can devise procedures which must neces- 
sarily converge at least to a local minimum, and 
which can be accelerated by the use of more sophis- 
ticated descent methods such as conjugate gradient 
and quasi-Newton algorithms. There is the possi- 
bility of more than one local minimum, but in any 
case the method will lead to an improvement over 
the original design. Furthermore, unlike the tradi- 
tional inverse algorithms, the cost function can be 
thought of as any measure of performance provided 
the adjoint system can be correctly formulated. 

In the present method the steps to obtain (1 - 4) 
are applied to the governing differential equations 
rather than applying them to the corresponding dis- 
crete system. This latter approach is now gaining fa- 
vor in the work of Newman and Tayor et al. [11, 8]. 
The two methods can be very similar, depending 
upon the discretization of (3). The current method 
has the advantage that the discretization and iter- 
ation scheme used to solve the flowfield system can 
also be applied directly to the adjoint system (4). In 
Jameson’s previous application of control theory to 
optimal aerodynamic shape design a successful nu- 
merical implementation was demonstrated using a 
conformal transformation from a unit circle to gen- 
erate the profile, so that the mapping function be- 
comes the control. In this work, an alternative ap- 
proach using a general coordinate transformation is 
adopted, and the equations are discretized by a finite 
volume method. This is intended to be a precursor 
for the three-dimensional problem, where conformal 
mapping is less suitable since it could be used only 
to provide independent transformations in separate 
planes. 


Design for potential flow using 
finite volume discretization 

Consider the case of two-dimensional compressible 
inviscid flow. In the absence of shock waves, an ini- 
tially irrotational flow will remain irrotational, and 
we can assume that the velocity vector q is the gra- 
dient of a potential <j>. In the presence of weak shock 
waves this remains a fairly good approximation. 
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Let p, p, c, and M be the pressure, density, speed- 
of-sound, and Mach number q/c. The potential flow 
equation is then 


V- {pVi) = 0, 

where the density is given by 
7-1 


P={l + L - r Ml (l-? 2 ) 


while 


P = 


c 2 = I2. 


(5) 


( 6 ) 


(7) 


iMl ' p 

Here M 0 0 is the Mach number in the free stream, 
and the equations have been non-dimensionalized so 
that p and q have the value unity in the far field. 
The potential flow equation can be written 

±(pu) + ±fa) = 0 inD. 

where u and t; represent the cartesian velocity com- 
ponents. The coordinate transformations may be 
defined 


u 


<t>x 


MU 
1 

in 

dx 


n 

V 


<8 y 


SL 
. a y 

in. 

dy _ 


. ^ . 


= H r ~ 


h 

4><, 


( 8 ) 


and also 




dx 

di 
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‘ <t>x ’ 

= h t 

<Px 



dx 



4>y 


_ <f> y 

L J 


„ 

dr 7 _ 






where x and y represent the physical plane, and Z 
and rj represent the computational plane. By defin- 
ing the Jacobian 

_ dx dy dx dy 

dZ di) drf dZ ’ 

we can write 

JL(phU) + -^(phV) = 0 in D. (9) 
Here, U and V represent the contravariant velocities 


dy dx 

dy drj 

dy dx 

di di 
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= H~ l H 2 
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<i>T, 


= J~ l 
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tv 


with 


Thus, 


J~ l = 


-Au 

-Al2 


An 

A 22 


U — An<f>£ -f Ai2<Pr, (10) 

V = Ai2<i>( + A22<Pr}- (11) 

Consider first the case in which the cost function 
is defined such as to achieve a target speed distribu- 


tion: 
I 


= \j c [q ~ qdf d3 = \j c [q - qd)i ($) 


( 12 ) 


where qd is the desired speed distribution and C is 
the airfoil surface. 

The design problem is now treated as a control 
problem where the control function is the airfoil 
shape, which is to be chosen to minimize I subject to 
the constraints defined by the flow equations (5-11). 
The first variation of this cost function is 

«- 


/ c <»-«>»» (5 

/('-«> I s < 13 > 


since on the wall 


d<f> 


dZ 

d<j> dZ 


qw ~ ds dZ ds ' 

In general we need to find how a modification to 
the airfoil geometry causes a variation 8<j>, as well 
a variation in the grid parameters 8A\\, 8 An, SA 2 2 , 
and 6h. Consider 

6U = <5 (.An) + A\\8<j>t + 8 (A 12 ) 1 fir) + An84>y 

8V = 8 (An) <8( + Au8<t>z -|- 8 (A 22 ) <t > >j + A 228 <Pn 


» = -5 




8$ 


~^~2 [ 6 A\\ 4 * + 28 An<?n<t>t + 8A22<t>~ r) \ ■ 


It follows that 8<j> satisfies 


U6 = 


pUSh 

+ Ph4>( 

+ph4n (l 8A 

( Aph^ SA 22 


8A n 

12 
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Finally the variation can be defined as 


pVSh 




+ ph<p n ( 1 - 


Vi* 

~2c 1 ' 

V*-, 


8 A 22 


+ ph<f>t ^1 — 8 An 
+ ph<p( (—• 8 An 
or, defining the right hand side with P and Q, 

L 8 <p = — -jr-Q ( 6 h, 5i4n, 8 An, 8 A 22 ) 

08 . 


with 


L = 


-—P(5h,8A xl ,8An,8A 22 ) (14) 
dr] 


ph (.An 

K[ + ph(An-^)f,\ 

r MAu-Wh 


d 

+ dr] 


[ +P h (^22 - p-) I; J 


• ( 15 ) 


If t/> is any periodic function vanishing in the far field, 
equation (14) can be multiplied by 4> and integrated 
oven- the domain. After integrating the right hand 
side by parts we arrive at 


J o *UW*l = J D % 


Q + j^P dtdr) 

drj 


■ J {i>ph[ 6 An<l>z + bAn^r,]} d£- (16) 


Now subtracting (16) from (13), 


61 = -J 

Jc 

6T. 


d(q- qd) 


8 <j>dZ 


+ 

+ 




d8 


{ * + #’*'* 


+ 


98 

■ j c \ (« - <ld ) 2 

I '(£)£« 

— [ xpL84>d£dTi 
Jo 

+ L 

• J {ii>ph {8An4>i + > 5 ^ 22 ^]} d£ : 

Then setting up the adjoint system we have 

£,0 = 0 in D, (17) 

with the boundary condition 

ph {Anrl>( + ^ 22 ^. 7 ) = — -q£ {q - q<t) • ( 18 ) 

After appling the second form of Green’s theorem to 
(17) we get 

f xpLS<j)dS = f {i>ph[8Au4>t + SAntq]} d£ 

Jd Jc 

+ J [8<pph ( 6 An^t + 5 ^ 22 ^ 17 ]} d£. 





For this finite volume formulation no general ana- 
lytic grid transformation is available. Furthermore, 
the variation with respect to the grid quantities is 
now spread into 6A U , 6 An, 8 A 22 . and 8 h instead of 
just the modulus of the transformation as was the 
case for conformal mapping. Therefore, we adopted 
a more conventional method to construct 81. First, 
an independent basis space of perturbation functions 
6(t), i = 1,2, ...,n (n = number of design vari- 
ables) is chosen that allows for the needed freedom 
of the design space. Thus, the shape / now becomes 
/(6(i)) with 6(i ) being the control. Once 81 is ob- 
tained, any optimization procedure can be employed 
to minimize the cost with respect to the given basis 

6(i). 

If the flow is subsonic, this procedure should con- 
verge toward the desired speed distribution since 
the solution will remain smooth, and no unbounded 
derivatives will appear. If, however, the flow is tran- 
sonic, one must allow for the appearance of shock 
waves in the trial solutions, even if qj is smooth. In 
such instances q — qd is not differentiable. This dif- 
ficulty can be circumvented by a more sophisticated 
choice of the cost function. Consider the choice 



where Aj and A 2 are parameters, and the periodic 
function S(() satisfies the equation 


. c . d*S 

Alh - A 2^2 = ? ~ Id- 


( 20 ) 


Then, 


SI = / c (w«S + A,|i«s)df 

= / S^SS-X^SsJ d( = J SSqdS. 

Thus, S replaces q — qt in the previous formula and 
one modifies the boundary condition (18) to 


ph (An^s + A 22 <f>ij) 


-A(S) onC. 


( 21 ) 


5 


For the case where the cost function is drag, (12) 
is replaced by, 


6. Project SI into a feasible direction subject to 
any constraints to obtain 61. 


l = < 22 > 

The first variation of the cost function is now, 



where the boundary condition on ip, (18) or (20), is 
replaced with 

nV'f +^ 22 4>n) = (25) 


d 2 S dy . . 

^ S ~ X ^=P q Ts' (26) 

The entire procedure can be summarized for the cost 
function based on target speed distribution as fol- 
lows: 


1. Solve the flow equations (5-11) for <p, u, v, q, p, 
U, and V. 

2. Smooth the cost function if necessary by (20). 

3. Solve the adjoint equation (15 and 17) for i p 
subject to the boundary condition (18) or (21). 

4. For each i independently perturb the design 
variables, 6(i), and calculate the necessary 
metric variations (Mu, 6Au, SAm, Sh, and 
^(sf)) by recalculating the perturbed grid 
with automatic gTid generation. 

5. Directly evaluate 61 by equation (19). 


7. Feed 61 as the gradient with respect to b(i) to 
a quasi-Newton optimization procedure. 

8. Calculate the search direction from the quasi- 
Newton algorithm and perform a line search. 

9. Return to 1 if not converged. 

In practice the method resembles those used by 
Hicks ef al. [13] with the control theory replacing 
the brute force, finite difference based, gradient cal- 
culation. The current formulation has an advan- 
tage by requiring computational work proportional 
to 2+m flow solver evaluations (m being the number 
of calculations required per line search) per design 
cycle as opposed to 1 + m + n. Thus, unlike con- 
ventional design optimization programs, the current 
method’s computational cost does not hinge upon 
the number of design variables provided the grid re- 
generation is fast and automatic. The method also 
has the advantage of being quite general in that ar- 
bitrary choices for both the design variables and op- 
timization technique are admitted. Finally, unlike 
the conformal mapping based method this approach 
can be directly extended to three-dimensions. 

Implementation of the generalized 
potential flow design method 

The practical implementation of the generalized po- 
tential flow design method, as with the conformal 
potential method, relies heavily upon fast accurate 
solvers for both the state (<£) and co-state (ip) fields. 
Further, to improve the speed and realizability of 
the methods, a robust choice of the optimization 
algorithm must be made. Finally, appropriate de- 
sign variables must be choosen which allow suffi- 
cient freedom in realizable designs. In this work, 
Jameson’s FL042 full potential computer program 
and the QNMDIF (by Gill, Murray and Wright [2]) 
quasi-Newton optimization algorithm are employed. 

In FL042 the flow solution is obtained by a rapid 
multigrid alternating direction method. The orig- 
inal scheme is described in [5]. The scheme uses 
artificial dissipative terms to introduce upwind bias- 
ing which simulates the rotated difference scheme [4] 
while preserving the conservation form. The alter- 
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nating direction method is a generalization of con- 
ventional alternating direction methods in which the 
scalar parameters are replaced by upwind difference 
operators to produce a scheme which remains stable 
as the equations change type from elliptic to hyper- 
bolic in accordance with the flow becoming locally 
supersonic [5]. 

QNMDIF is an unconstrained quasi-Newton 
optimization algorithm that calculates updates 
to a Cholesky factored Hessian matrix by the 
DFP (Davidon-Fletcher-Powell) rank-two proce- 
dure. Hence, information about the curvature of the 
design space feeds in through the successive gradient 
calculations. 

Since the primary computational costs arise from 
not only the flow solution algorithm but also the 
adjoint solution algorithm, both need to be compu- 
tationally efficient. The adjoint equation has a form 
very similar to the flow equation. While it is linear in 
its dependent variable, it also changes type from el- 
liptic (in subsonic zones of the flow) to hyperbolic (in 
supersonic zones of the flow). Thus, it was possible 
to adapt exactly the same algorithm to solve both 
the adjoint and the flow equations, but with reverse 
biasing of the difference operators in the downwind 
direction for the adjoint equation, corresponding to 
its reversed direction of the zone of dependence. A 
multigrid method is used to accelerate the conver- 
gence of a generalized alterating direction scheme in 
a manner similar to the flow solver. 

Design variables are chosen with the following 
form, suggested by Hicks and Henne [3]: 


b (z) = sin 



i(z) S i'‘(L-:)e* ,,t , 

where tj and *2 control the center and thickness 
of the perturbation and z is the normalized chord 
length. 

When distributed over the entire chord on both 
upper and lower surfaces these analytic perturbation 
functions admit a large possible design space. They 
have the advantage of being space based functions, 
as opposed to frequency based functions, and thus 
they allow for local control of the design. They can 
be chosen such that symetry, thickness, or volume 
can be explicitly constrained. Further, particular 
choices of these variables will concentrate the design 


effort in regions where refinement is needed, while 
leaving the rest of the airfoil section virtually undis- 
turbed. The disadvantage of these functions is that 
they do not form a complete basis space, nor are 
they orthogonal. Thus, they do not garantee that a 
solution, for example, of the inverse problem for a 
realizable target pressure distribution will necessar- 
ily be attained. Here they are employed due to their 
ease of use and ability to produce a wide variation 
of shapes with a limited number of design variables. 


Numerical tests of the generalized 
potential flow design method 


Several test cases are presented for the generalized 
potential flow design algorithm based on the finite 
volume scheme. These test cases can be catagorized 
into three basic groups. First, are non-lifting cases, 
where a symetric target pressure distribution is spec- 
ified and the optimization is started from an arbi- 
trary symetric initial guess. 

The first non-lifting example shown in Figure 1, 
illustrates that for subsonic flow, M = 0.2 and a = 
0°, a given airfoil shape, in this case a NACA 64012, 
can be recovered by starting from an arbitrary shape 
and specifying the target pressure distribution. A 
close look at the final solution shows that a small 
discrepency is evident at the trailing edge. This 
may be associated with the lack of completness of 
our basis space. In the next example, see Figure 2, 
the design takes place at Moo = 0.8, a — 0°, where 
the initial NACA 0012 airfoil is driven towards the 
subsonic pressure distribution of the NACA 64021. 
In this case the target pressure distribution exceeds 
Cp* for Moo = 0-8- Therefore, the pressure distri- 
bution represents shock free transonic flow. Since, 
in general, such a pressure distribution may not be 
realizable, the program approaches the target with 
the nearest feasible pressure distribution. An exam- 
ination of Figure 2 demonstrates that a very weak 
shock in the designed pressure distribution replaces 
the smooth transition to subsonic flow seen in the 
target distribution. In the final example non-lifting 
case of Figure 3, an arbitrary pressure distribution 
which does contain a shock wave and is realizable, 
is used as the target. Here the computer program 
was able to obtain the corresponding airfoil geome- 
try along with the correct shock wave location with 
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a high degree of accuracy, as can be seen both in the 
pressure distribution and in the airfoils. 

The second group of test cases address the prob- 
lem of attaining a desired pressure distribution for 
lifting airfoils. The most convienient method of 
obtaining such solutions with the present design 
method is to determine the lift coefficient associated 
with the target pressure distribution, and match this 
lift with the initial airfoil. The design progresses 
with the flow solver and the adjoint system being 
driven by constant circulation instead of fixed an- 
gle of attack. The first example using this tech- 
nique, shown in Figure 4, drives the NACA 0012 
airfoil toward the target pressure distribution for 
the NACA 64A410 airfoil at M „ = 0.735, a = 0°, 
and C( = 0.75. This case requires a shift in the 
shock location and a significant change in the pro- 
file shape such that the target pressure distribution 
is obtained. The final solution almost exactly recov- 
ers the pressure distribution and the airfoil shape. 
In the next example, Figure 5, the NACA 0012 air- 
foil is again used as the starting condition to obtain 
the pressure distribution of the GAW72 airfoil oper- 
ating at Mo, = 0.7, a - -2°, and C/ = 0.57. This 
case is difficult since the target airfoil has a cusped 
trailing edge while the initial airfoil has a finite trail- 
ing edge. As was seen in some of the non-lifting 
cases, there are small discrepencies evident near the 
trailing edge that may be due to the incomplete ba- 
sis of the chosen design variables. The difference 
in the profiles between the final design and actual 
GAW72 is partly due to the fact that the GAW72 
cordinates place the trailing edge at a non- zero y or- 
dinate while the NACA 0012 places the trailing edge 
at y = 0. Also, the designed airfoil is subject to an 
arbritrary rotation since the angle of attack is free 
during optimization. The last test case in which the 
design program is run in inverse mode involves driv- 
ing the NACA 0012 airfoil at Moo = 0.75 to obtain 
the target pressure distribution of the RAE airfoil 
at the same Mach number, a = 1.0°, and Cj = 0.80. 
Due to the steep favorable pressure grardient at the 
leading edge upper surface and the strong shock ex- 
hibited (see Figure 6) by the RAE airfoil at these 
conditions this case represents quite a difficult test 
for the program. In the observed results, discrepen- 
cies are evident between the target and the final de- 
signed pressure distributions in both the leading and 


trailing edge regions. A comparison of the profiles 
reveals that the floating a in the design process has 
resulted in a rotation between the target and the 
final design. Looking at the final angle of attack 
of 0.74° which is slightly off the a = 1.0° of the 
target reveals that the rotation is compensated in 
large degree by the floating a. However, both air- 
foils have leading and trailing edges at 0.0; hence, 
the difference between the airfoils is not a simple 
rotation. This result, combined with an incomplete 
basis space, may account for the observed differences 
in pressure distributions. 

The last group of results introduces drag as the 
cost function. Again the design process is carried 
out in the fixed lift mode. In Figure 7, the first 
drag minimization example, a NACA 0012 is again 
used as a starting airfoil. The design takes place at 
M 0 o = 0.75 and Ci = 0.50 where a strong shock 
causes considerable wave drag in the initial airfoil. 
To make the problem interesting, the optimization 
is carried out such that symetry of the design is 
perserved. The final design is a symetric airfoil with 
an increased maximum thickness that operates at 
the same lift coefficient, but has a reduction in drag 
from Cd = 0.0127 to Cd = 0.0016. In the final 
test case (see Figure 6) the camber distribution is 
optimized insted of thickness distribution. The de- 
sign starts from a NACA 64A410 airfoil operating at 
Moo = 0.75, and Ci = 0.60 which displays 42 counts 
of drag according to the potential flow calculation. 
By allowing only changes to the camber distribution, 
a final airfoil is produced which maintains C/ = 0.60 
but does so with only 4 counts of drag. 

Conclusions and recommendations 

We have developed a control theory based airfoil 
design method for a two-dimensional finite volume 
discretization of the potential flow equation. The 
method represents an extension of Jameson’s pre- 
vious work on the design of a conformally mapped 
airfoil. The new method is both efficient and robust, 
combining the versatility of numerical optimization 
methods with the efficiency of inverse design meth- 
ods. The motivating factor behind this work is its 
direct extendability to three-dimensions. The ulti- 
mate goal of this effort is to create practical aero- 
dynamic shape design methods for complete aircraft 
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' configurations. 
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la: Initial Condition 
C, = 0.0000, C d = 0.0001 


lb: 7 Design Iterations 
C, = 0.0000, C d = 0.0000 


Figure 1: Subsonic Non-Lifting Design Case, M = 0.2, a = 0°. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p : NACA 64012, M = 0.2. 

Inverse Design 



2a: Initial Condition 
C, = 0.0000, C d = 0.0063 



C 



2b: 7 Design Iterations 
C, = 0.0000, C d = 0.0003 


Figure 2: Transonic Non-Lifting Design Case, M = 0.8 a = 0°. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p : NACA 64021, M = 0.2. 

Inverse Design 
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3a: Initial Condition 
Ci = 0.0000, C d = 0.0063 



3b: 8 Design Iterations 
C, = 0.0000, Cd = 0.0015 


Figure 3: Transonic Non-Lifting Design Case, M = 0.8, a = 0°. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C ? : NACA 64X, M = 0.8. 

Inverse Design 



Figure 4: Transonic Lifting Design Case, M = 0.735 Fixed Lift. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p : NACA 64A410, M = 0.735, C, = 0.73. 
Inverse Design 
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5a: Initial Condition 
C, = 0.5492, C d = 0.0047, a = 2.709° 



C, = 0.5496, C d = 0.0045, a = -1.508° 


Figure 5: Transonic Lifting Design Case, M = 0.70, Fixed Lift. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p : GAW72, M = 0.70. 

Inverse Design 



Figure 6: Transonic Lifting Design Case, M = 0.75 Fixed Lift. 
— , x Initial Airfoil: NACA 0012. 

, + Target C p : RAE, M = 0.75. 

Inverse Design 
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Figure 7: Transonic Lifting Design Case, M — 0.75, Fixed Lift. 
— , x Initial Airfoil: NACA 0012. 

Symetric Drag Minimization. 



8a: Initial Condition 
Ci = 0.5964, Ci = 0.0042, a = -0.464° 



8b: 2 Design Iterations 
C, = 0.5966, Ci = 0.000*, a = 0.175° 



Figure 8: Transonic Lifting Design Case, M = 0.735 Fixed Lift. 
— , x Initial Airfoil: NACA 64A410. 

Camber Only Drag Minimization. 
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